/*
 * NewtonSolver.cpp
 *
 *  Created on: 20 Oct 2011
 *      Author: allan
 */

#include "NewtonSolver.h"

// GeoReact includes
#include <Utils/Assert.h>

NewtonSolver::NewtonSolver() :
decomposition(FullPivotingLU), tolerance(1.0E-6), maxiterations(100)
{}

NewtonSolver::NewtonSolver(unsigned dimension) :
decomposition(FullPivotingLU), tolerance(1.0E-6), maxiterations(100)
{
	Initialiaze(dimension);
}

void NewtonSolver::Initialiaze(unsigned dimension)
{
	F = VectorXd::Zero(dimension);
	J = MatrixXd::Zero(dimension, dimension);
}

void NewtonSolver::SetDecompositionAlgorithm(Decomposition decomposition)
{
	this->decomposition = decomposition;
}

void NewtonSolver::SetConvergenceOptions(double tolerance, uint maxiterations)
{
	this->tolerance = tolerance;
	this->maxiterations = maxiterations;
}
	
bool NewtonSolver::Iterate(NewtonProblem& problem, VectorXd& x)
{
	problem.Function(x, F);
	problem.Jacobian(x, J);
	
	VectorXd dx;
	
	switch(decomposition)
	{
	case PartialPivotingLU:
		     partial_lu.compute(J);
		dx = partial_lu.solve(-F);
		break;
	default:
		     full_lu.compute(J);
		dx = full_lu.solve(-F);
		break;
	}
	
	x += dx;
	
	Assert((J * dx).isApprox(-F), "The linear system could not be solved or does not have a solution because of singular matrix.");
	
	return (F.norm() < tolerance);
}

ConvergenceDiagnostics NewtonSolver::Solve(NewtonProblem& problem, VectorXd& x)
{
	ConvergenceDiagnostics diagnostics = {0, false};
	
	while(!diagnostics.converged && diagnostics.iterations++ < maxiterations) 
		diagnostics.converged = Iterate(problem, x);
	
	return diagnostics;
}
